import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
import healpy as hp
import pixell as px
from pixell import reproject
from classy import Class

pi = np.pi

plt.style.use(['mysize_meeting'])
fig = plt.figure()
fw = fig.get_figwidth()
fh = fig.get_figheight()

<Figure size 720x540 with 0 Axes>

cmap = plt.cm.RdBu

c__256 = 'C0'
c__512 = 'C1'
c_1024 = 'C2'
c_proj = 'C2'
c_fgen = 'C3'

 

 

Size and resolution

The resolution of a full-sky map is governed by the parameter nside. This is also where I specify the image size used throughout this notebook.

res2048_arcmin = hp.nside2resol(2048, arcmin=True)
res1024_arcmin = hp.nside2resol(1024, arcmin=True)
res512_arcmin = hp.nside2resol(512, arcmin=True)
res256_arcmin = hp.nside2resol(256, arcmin=True)
print(res2048_arcmin)
print(res1024_arcmin)
print(res512_arcmin)
print(res256_arcmin)

1.717743205908703
3.435486411817406
6.870972823634812
13.741945647269624

imgsize_deg = 100  # in degrees
imgsize_arcmin = imgsize_deg * 60  # in arcmin
imgsize_rad = imgsize_deg / 180 * pi  # in radians
pixsize_256  = int(np.ceil(imgsize_arcmin / res256_arcmin))
pixsize_512  = int(np.ceil(imgsize_arcmin / res512_arcmin))
pixsize_1024 = int(np.ceil(imgsize_arcmin / res1024_arcmin))
imgsize_arcmin = pixsize_1024 * res1024_arcmin
print("img: %.2f deg = %.f arcmin" % (imgsize_arcmin/60, imgsize_arcmin))
print("pix: %d pixels" % pixsize_1024)
print("res: %.3f deg = %.2f arcmin" % (res1024_arcmin/60, res1024_arcmin))

img: 100.03 deg = 6002 arcmin
pix: 1747 pixels
res: 0.057 deg = 3.44 arcmin

The theoretical input power spectrum

from the CMB

try:
    cosmo.struct_cleanup()
except:
    pass
cosmo = Class()
cosmo.set({
    'output': 'tCl lCl', 
    'modes': 's',
    'omega_b': 0.022383, 
    'omega_cdm': 0.12011, 
    '100*theta_s': 1.041085, 
    'tau_reio': 0.0543, 
    'A_s': np.exp(3.0448)*1e-10, 
    'n_s': 0.96605, 
    'm_ncdm': 0.06, 
    'N_ncdm': 1, 
    'N_ur': 2.0328, 
    'lensing': 'yes', 
    'l_max_scalars': 20000
})

True

%%time
cosmo.compute()

CPU times: user 1min 30s, sys: 3.28 s, total: 1min 33s
Wall time: 7.33 s

ell = cosmo.lensed_cl()['ell']
lmax = ell[-1]
norm = ell * (ell+1) / (2*np.pi)
cl_cmb = cosmo.lensed_cl()['tt'] * cosmo.T_cmb()**2 / 1e-6**2

from custom power spectrum

cl_pwl = ell**2 / (1e9 + ell**4) * 2e8
cl_pwl[1:] /= norm[1:]

Plot

def power_axes(ax1, ax2):
    ax1.set_xlabel("$\ell$")
    ax2.set_xlabel("$\ell$")
    ax1.set_ylabel("$C_\ell\ \ell (\ell+1) / 2\pi\ /\ \mu\mathrm{K^2}$")
    ax2.set_ylabel("$C_\ell\ \ell (\ell+1) / 2\pi\ /\ \mu\mathrm{K^2}$")
    ax1.set_xlim(-100, 3000)
    ax2.set_xlim(1, 10000)

fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(ell, cl_cmb*norm, c='k',   label="CMB power spectrum (from CLASS with P18 best-fit parameters)")
ax1.plot(ell, cl_pwl*norm, c='0.5', label="custom power-laws power spectrum")
ax1.legend(handlelength=0, labelcolor='linecolor')

ax2 = plt.subplot(212)
ax2.loglog(ell, cl_cmb*norm, c='k'  )
ax2.loglog(ell, cl_pwl*norm, c='0.5')
ax2.set_ylim((cl_pwl*norm)[2]/10, (cl_pwl*norm).max()*4)

power_axes(ax1, ax2)

 

From angular power spectrum to map

cl = cl_cmbcl = cl_peak

%%time
alm_cmb = hp.synalm(cl_cmb, lmax=lmax)
alm_pwl = hp.synalm(cl_pwl, lmax=lmax)

CPU times: user 30.5 s, sys: 1.8 s, total: 32.3 s
Wall time: 33.8 s

%%time
map_cmb_256  = hp.alm2map(alm_cmb, nside=256,  lmax=lmax)
map_cmb_512  = hp.alm2map(alm_cmb, nside=512,  lmax=lmax)
map_cmb_1024 = hp.alm2map(alm_cmb, nside=1024, lmax=lmax)
map_pwl_256  = hp.alm2map(alm_pwl, nside=256,  lmax=lmax)
map_pwl_512  = hp.alm2map(alm_pwl, nside=512,  lmax=lmax)
map_pwl_1024 = hp.alm2map(alm_pwl, nside=1024, lmax=lmax)

CPU times: user 15min 14s, sys: 705 ms, total: 15min 14s
Wall time: 1min 1s

vmax = max(abs(map_cmb_1024.min()), abs(map_cmb_1024.max()))
hp.mollview(map_cmb_1024, min=-vmax, max=vmax, cmap=cmap, title="CMB map")

vmax = max(abs(map_pwl_1024.min()), abs(map_pwl_1024.max()))
hp.mollview(map_pwl_1024, min=-vmax, max=vmax, cmap=cmap, title="custom power-law map")

 

From map to angular power spectrum

%%time
cl_cmb_out_256  = hp.anafast(map_cmb_256,  lmax=lmax, use_pixel_weights=True)
cl_cmb_out_512  = hp.anafast(map_cmb_512,  lmax=lmax, use_pixel_weights=True)
cl_cmb_out_1024 = hp.anafast(map_cmb_1024, lmax=lmax, use_pixel_weights=True)
cl_pwl_out_256  = hp.anafast(map_pwl_256,  lmax=lmax, use_pixel_weights=True)
cl_pwl_out_512  = hp.anafast(map_pwl_512,  lmax=lmax, use_pixel_weights=True)
cl_pwl_out_1024 = hp.anafast(map_pwl_1024, lmax=lmax, use_pixel_weights=True)

CPU times: user 15min 7s, sys: 3.72 s, total: 15min 11s
Wall time: 1min 13s

fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(ell, cl_cmb_out_256  * norm, c=c__256, label="from map with nside=256")
ax1.plot(ell, cl_cmb_out_512  * norm, c=c__512, label="from map with nside=512")
ax1.plot(ell, cl_cmb_out_1024 * norm, c=c_1024, label="from map with nside=1024")
ax1.plot(ell, cl_cmb * norm, c='k', label="CMB input Cl")
ax1.set_ylim(0, max(cl_cmb*norm)*1.2)
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor')

ax2 = plt.subplot(212)
ax2.loglog(ell, cl_cmb_out_256  * norm, c=c__256)
ax2.loglog(ell, cl_cmb_out_512  * norm, c=c__512)
ax2.loglog(ell, cl_cmb_out_1024 * norm, c=c_1024)
ax2.loglog(ell, cl_cmb * norm, c='k')
ax2.set_ylim((cl_cmb*norm)[2]*1e-2, max(cl_cmb*norm)*4)

power_axes(ax1, ax2)

fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(ell, cl_pwl_out_256  * norm, c=c__256, label="from map with nside=256")
ax1.plot(ell, cl_pwl_out_512  * norm, c=c__512, label="from map with nside=512")
ax1.plot(ell, cl_pwl_out_1024 * norm, c=c_1024, label="from map with nside=1024")
ax1.plot(ell, cl_pwl * norm, c='k', label="custom power-law input Cl")
ax1.set_ylim(0, max(cl_pwl*norm)*1.2)
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor')

ax2 = plt.subplot(212)
ax2.loglog(ell, cl_pwl_out_256  * norm, c=c__256)
ax2.loglog(ell, cl_pwl_out_512  * norm, c=c__512)
ax2.loglog(ell, cl_pwl_out_1024 * norm, c=c_1024)
ax2.loglog(ell, cl_pwl * norm, c='k')
ax2.set_ylim((cl_pwl*norm)[2]*1e-2, max(cl_pwl*norm)*4)

power_axes(ax1, ax2)

Fourier power spectrum

def fourier_power_spectrum(img_r, img_len, bin_k):
    """
    Compute the Fourier power spectrum given a 2-dimensional real-space map.
    
    Parameters
    ----------
        img_r: arraylike
            2-dimensional real-space map 
            
        img_len: float
            image size 
            (units define output units of wavenumber and power spectrum)
            
        bin_k: int
            number of bins in Fourier space
            
    Returns
    -------
        k: np.ndarray
            wavenumber in units of inverse [img_len]
            
        pk: np.ndarray
            power spectrum in units of [img_len]**2
    """
    assert np.ndim(img_r) == 2, 'img_r is not 2D'
    img_pix = np.size(img_r, axis=0)  # pixel number per axis
    img_dim = np.ndim(img_r)  # image dimension

    # This first 'paragraph' is to create masks of indices corresponding to 
    # one Fourier bin each.
    freq = np.fft.fftfreq(n=img_pix, d=img_len/img_pix) * 2*np.pi
    rfreq = np.fft.rfftfreq(n=img_pix, d=img_len/img_pix) * 2*np.pi
    kx, ky = np.meshgrid(freq, rfreq, indexing='ij')
    k_abs = np.sqrt(kx**2. + ky**2.)
    # The following complicated line is actually only creating a 1D array 
    # spanning k-space logarithmically from minimum k_abs to maximum. To 
    # start slightly below the minimum and finish slightly above the maximum 
    # I use ceil and floor. To ceil and floor not to the next integer but to 
    # the next 15th digit, I multiply by 1e15 before flooring and divide by 
    # the same amount afterwards. Since the ceiled/floored value is actually 
    # the exponent used for the logspace, going to the next integer would be 
    # too much, which is why I do this ceil/floor to digit procedure.
    k_log = np.logspace(np.floor(np.log10(np.min(k_abs[1:]))*1.e15)/1.e15, np.ceil(np.log10(np.max(k_abs[1:]))*1.e15)/1.e15, bin_k)

    img_k = np.fft.rfftn(np.fft.fftshift(img_r)) * (img_len/img_pix)**img_dim
    pk = np.empty(np.size(k_log)-1)
    for i in range(np.size(k_log)-1):
        mask = (k_abs >= k_log[i]) & (k_abs < k_log[i+1])
        pk[i] = np.mean(np.abs(img_k[mask])**2.) / img_len**img_dim

    k = (k_log[1:] + k_log[:-1]) / 2

    return k, pk

 

Projection to flat map

Using the pixell package as mentioned in the ACT example notebooks to project the full-sky map to a flat cut-out.

from pixell import enmap, utils

box = np.array([[-imgsize_arcmin//2, imgsize_arcmin//2], [imgsize_arcmin//2, -imgsize_arcmin//2]]) * utils.arcmin
shape256, wcs256 = enmap.geometry(pos=box, res=res256_arcmin*utils.arcmin, proj='car')
shape512, wcs512 = enmap.geometry(pos=box, res=res512_arcmin*utils.arcmin, proj='car')
shape1024, wcs1024 = enmap.geometry(pos=box, res=res1024_arcmin*utils.arcmin, proj='car')
imap256  = enmap.zeros((3,) + shape256,  wcs=wcs256)
imap512  = enmap.zeros((3,) + shape512,  wcs=wcs512)
imap1024 = enmap.zeros((3,) + shape1024, wcs=wcs1024)

imap1024.shape

(3, 1747, 1747)

imap1024.wcs

car:{cdelt:[-0.05726,0.05726],crval:[-0.008333,0],crpix:[874.1,874]}

%%time
fmap_cmb_256  = reproject.enmap_from_healpix(map_cmb_256,  shape=(pixsize_256,  pixsize_256),  wcs=imap256.wcs)
fmap_cmb_512  = reproject.enmap_from_healpix(map_cmb_512,  shape=(pixsize_512,  pixsize_512),  wcs=imap512.wcs)
fmap_cmb_1024 = reproject.enmap_from_healpix(map_cmb_1024, shape=(pixsize_1024, pixsize_1024), wcs=imap1024.wcs)
fmap_pwl_256  = reproject.enmap_from_healpix(map_pwl_256,  shape=(pixsize_256,  pixsize_256),  wcs=imap256.wcs)
fmap_pwl_512  = reproject.enmap_from_healpix(map_pwl_512,  shape=(pixsize_512,  pixsize_512),  wcs=imap512.wcs)
fmap_pwl_1024 = reproject.enmap_from_healpix(map_pwl_1024, shape=(pixsize_1024, pixsize_1024), wcs=imap1024.wcs)

Preparing SHT
T -> alm
float64 complex128
Computing pixel positions
Computing rotated positions
Projecting
Preparing SHT
T -> alm
float64 complex128
Computing pixel positions
Computing rotated positions
Projecting
Preparing SHT
T -> alm
float64 complex128
Computing pixel positions
Computing rotated positions
Projecting
Preparing SHT
T -> alm
float64 complex128
Computing pixel positions
Computing rotated positions
Projecting
Preparing SHT
T -> alm
float64 complex128
Computing pixel positions
Computing rotated positions
Projecting
Preparing SHT
T -> alm
float64 complex128
Computing pixel positions
Computing rotated positions
Projecting
CPU times: user 57.5 s, sys: 980 ms, total: 58.5 s
Wall time: 13.7 s

vmax = max(abs(fmap_cmb_1024[0].min()), abs(fmap_cmb_1024[0].max()))
plt.title("flat projection from CMB map\n %.f x %.f deg$^2$ \n %d x %d pixels" 
          % (imgsize_arcmin/60, imgsize_arcmin/60, pixsize_1024, pixsize_1024))
plt.imshow(fmap_cmb_1024[0], cmap=cmap, vmin=-vmax, vmax=vmax)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fa20efbaa90>

vmax = max(abs(fmap_pwl_1024[0].min()), abs(fmap_pwl_1024[0].max()))
plt.title("flat projection from custom power-law map\n %.f x %.f deg$^2$ \n %d x %d pixels" 
          % (imgsize_arcmin/60, imgsize_arcmin/60, pixsize_1024, pixsize_1024))
plt.imshow(fmap_pwl_1024[0], cmap=cmap, vmin=-vmax, vmax=vmax)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fa20f0c1040>

Fourier spectra

fmap_cmb_256_k,  fmap_cmb_256_pk  = fourier_power_spectrum(fmap_cmb_256[0],  img_len=imgsize_arcmin, bin_k=256)
fmap_cmb_512_k,  fmap_cmb_512_pk  = fourier_power_spectrum(fmap_cmb_512[0],  img_len=imgsize_arcmin, bin_k=256)
fmap_cmb_1024_k, fmap_cmb_1024_pk = fourier_power_spectrum(fmap_cmb_1024[0], img_len=imgsize_arcmin, bin_k=256)
fmap_pwl_256_k,  fmap_pwl_256_pk  = fourier_power_spectrum(fmap_pwl_256[0],  img_len=imgsize_arcmin, bin_k=256)
fmap_pwl_512_k,  fmap_pwl_512_pk  = fourier_power_spectrum(fmap_pwl_512[0],  img_len=imgsize_arcmin, bin_k=256)
fmap_pwl_1024_k, fmap_pwl_1024_pk = fourier_power_spectrum(fmap_pwl_1024[0], img_len=imgsize_arcmin, bin_k=256)

/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice.
  return _methods._mean(a, axis=axis, dtype=dtype,
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice.
  return _methods._mean(a, axis=axis, dtype=dtype,
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice.
  return _methods._mean(a, axis=axis, dtype=dtype,
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice.
  return _methods._mean(a, axis=axis, dtype=dtype,
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice.
  return _methods._mean(a, axis=axis, dtype=dtype,
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice.
  return _methods._mean(a, axis=axis, dtype=dtype,
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)

def power_axes_k(ax1, ax2):
    ax1.set_xlabel("$k$")
    ax2.set_xlabel("$k$")
    ax1.set_ylabel("$P(k)\ k^2 / (2\pi)^2$")
    ax2.set_ylabel("$P(k)\ k^2 / (2\pi)^2$")

fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(fmap_cmb_1024_k, fmap_cmb_1024_k**2/(2*pi)**2 * fmap_cmb_1024_pk, '-o', ms=2, c=c_1024, label="nside=1024")
ax1.plot(fmap_cmb_512_k, fmap_cmb_512_k**2/(2*pi)**2 * fmap_cmb_512_pk, '-o', ms=2, c=c__512, label="nside=512")
ax1.plot(fmap_cmb_256_k, fmap_cmb_256_k**2/(2*pi)**2 * fmap_cmb_256_pk, '-o', ms=2, c=c__256, label="nside=256")
ax1.set_ylim(min(np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10, 
             np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10))
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor', title="Fourier power spectrum \nfrom the $\\bf{projection}$ to a flat map \nof the CMB power spectrum")

ax2 = plt.subplot(212)
ax2.loglog(fmap_cmb_1024_k, fmap_cmb_1024_k**2/(2*pi)**2 * fmap_cmb_1024_pk, '-o', ms=2, c=c_1024)
ax2.loglog(fmap_cmb_512_k, fmap_cmb_512_k**2/(2*pi)**2 * fmap_cmb_512_pk, '-o', ms=2, c=c__512)
ax2.loglog(fmap_cmb_256_k, fmap_cmb_256_k**2/(2*pi)**2 * fmap_cmb_256_pk, '-o', ms=2, c=c__256)
ax2.set_ylim(min(np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10, 
             np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10))

power_axes_k(ax1, ax2)

fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(fmap_pwl_1024_k, fmap_pwl_1024_k**2/(2*pi)**2 * fmap_pwl_1024_pk, '-o', ms=2, c=c_1024, label="nside=1024")
ax1.plot(fmap_pwl_512_k, fmap_pwl_512_k**2/(2*pi)**2 * fmap_pwl_512_pk, '-o', ms=2, c=c__512, label="nside=512")
ax1.plot(fmap_pwl_256_k, fmap_pwl_256_k**2/(2*pi)**2 * fmap_pwl_256_pk, '-o', ms=2, c=c__256, label="nside=256")
ax1.set_ylim(min(np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10, 
             np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10))
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor', title="Fourier power spectrum \nfrom the $\\bf{projection}$ to a flat map \nof the custom power spectrum")

ax2 = plt.subplot(212)
ax2.loglog(fmap_pwl_1024_k, fmap_pwl_1024_k**2/(2*pi)**2 * fmap_pwl_1024_pk, '-o', ms=2, c=c_1024)
ax2.loglog(fmap_pwl_512_k, fmap_pwl_512_k**2/(2*pi)**2 * fmap_pwl_512_pk, '-o', ms=2, c=c__512)
ax2.loglog(fmap_pwl_256_k, fmap_pwl_256_k**2/(2*pi)**2 * fmap_pwl_256_pk, '-o', ms=2, c=c__256)
ax2.set_ylim(min(np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10, 
             np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10))

power_axes_k(ax1, ax2)

 

 

Generate Gaussian fluctuations on a flat map

class GaussianFluctuationMap:
    """
    Generator class for anisotropy maps from of Gaussian fluctuations for given angular power spectra.
    """
    def __init__(self, ell, cl, len_rad, res_rad):
        """
        Parameters
        ----------
            ell: list or np.ndarray
                multipole moment (needs to have the same shape as cl)
                
            cl: list or np.ndarray
                angular power spectrum (needs to have the same shape as ell)
                
            len_rad: float
                desired image size in radians
                
            res_rad: float
                desired image resolution in radians
        """
        self.len_rad = len_rad                # size of map in radians on the sky
        self.res_rad = res_rad                # resolution of the map in radians
        self.len_pix = self.rad2pix(len_rad)  # size of map in pixels

        self.ell = ell[2:]
        self.cl = cl[2:]
        self.l2cl = interp1d(ell, cl)

    def rad2pix(self, rad):
        """
        Convert radians in sky to pixels in an image at a given resolution.
        
        Parameters
        ----------
            rad: float or np.ndarray
            
        Returns
        -------
            pix: float or np.ndarray
                number of pixels corresponding to the given angel `rad` at the given resolution
        """
        pix = rad / self.res_rad
        if type(pix) == float:
            return int(pix)
        else:
            return pix.astype(int)

    @staticmethod
    def k2l_rad(k_irad):
        """
        Convert wavenumber k to multipolemoment l.
        
        Parameters
        ----------
            k_irad: float or np.ndarray
                wavenumber k in inverse radians
                
        Returns
        -------
            multipolemoment `ell` (dimensionless)
        """
        return 2 * pi * k_irad

    def cl2map(self):
        """
        Generate a random realization from the given specifications.
        """
        k_rfreq = np.fft.rfftfreq(self.len_pix, d=self.res_rad)
        k_freq = np.fft.fftfreq(self.len_pix, d=self.res_rad)
        kx, ky = np.meshgrid(k_freq, k_rfreq, indexing='ij')
        k_abs = np.sqrt(kx**2 + ky**2)
        cl = self.l2cl(self.k2l_rad(k_abs))
        # create the random map in k-space:
        k_img = np.sqrt(cl * self.len_rad**2/2) * (np.random.randn(*cl.shape) + 
                                              1j * np.random.randn(*cl.shape))
        # convert the map to real space (using an inverse Fourier transform):
        r_img = np.fft.fftshift(np.fft.irfft2(k_img)) / self.res_rad**2
        # devision by res to account for dx dy in FT, **2 for 2D
        return r_img

gm_cmb_256  = GaussianFluctuationMap(ell=ell, 
                                     cl=cl_cmb, 
                                     len_rad=imgsize_arcmin/60/180*pi, 
                                     res_rad=res256_arcmin/60/180*pi)
gm_cmb_512  = GaussianFluctuationMap(ell=ell, 
                                     cl=cl_cmb, 
                                     len_rad=imgsize_arcmin/60/180*pi, 
                                     res_rad=res512_arcmin/60/180*pi)
gm_cmb_1024 = GaussianFluctuationMap(ell=ell, 
                                     cl=cl_cmb, 
                                     len_rad=imgsize_arcmin/60/180*pi, 
                                     res_rad=res1024_arcmin/60/180*pi)
gm_pwl_256  = GaussianFluctuationMap(ell=ell, 
                                     cl=cl_pwl, 
                                     len_rad=imgsize_arcmin/60/180*pi, 
                                     res_rad=res256_arcmin/60/180*pi)
gm_pwl_512  = GaussianFluctuationMap(ell=ell, 
                                     cl=cl_pwl, 
                                     len_rad=imgsize_arcmin/60/180*pi, 
                                     res_rad=res512_arcmin/60/180*pi)
gm_pwl_1024 = GaussianFluctuationMap(ell=ell, 
                                     cl=cl_pwl, 
                                     len_rad=imgsize_arcmin/60/180*pi, 
                                     res_rad=res1024_arcmin/60/180*pi)

Map generation

gmap_cmb_256  = gm_cmb_256.cl2map()
gmap_cmb_512  = gm_cmb_512.cl2map()
gmap_cmb_1024 = gm_cmb_1024.cl2map()
gmap_pwl_256  = gm_pwl_256.cl2map()
gmap_pwl_512  = gm_pwl_512.cl2map()
gmap_pwl_1024 = gm_pwl_1024.cl2map()

Fourier spectra

gmap_cmb_256_k,  gmap_cmb_256_pk  = fourier_power_spectrum(gmap_cmb_256,  img_len=imgsize_arcmin, bin_k=512)
gmap_cmb_512_k,  gmap_cmb_512_pk  = fourier_power_spectrum(gmap_cmb_512,  img_len=imgsize_arcmin, bin_k=512)
gmap_cmb_1024_k, gmap_cmb_1024_pk = fourier_power_spectrum(gmap_cmb_1024, img_len=imgsize_arcmin, bin_k=512)
gmap_pwl_256_k,  gmap_pwl_256_pk  = fourier_power_spectrum(gmap_pwl_256,  img_len=imgsize_arcmin, bin_k=512)
gmap_pwl_512_k,  gmap_pwl_512_pk  = fourier_power_spectrum(gmap_pwl_512,  img_len=imgsize_arcmin, bin_k=512)
gmap_pwl_1024_k, gmap_pwl_1024_pk = fourier_power_spectrum(gmap_pwl_1024, img_len=imgsize_arcmin, bin_k=512)

/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice.
  return _methods._mean(a, axis=axis, dtype=dtype,
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)

fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(gmap_cmb_1024_k, gmap_cmb_1024_k**2/(2*pi)**2 * gmap_cmb_1024_pk, '-o', ms=2, c=c_1024, label="nside=1024")
ax1.plot(gmap_cmb_512_k, gmap_cmb_512_k**2/(2*pi)**2 * gmap_cmb_512_pk, '-o', ms=2, c=c__512, label="nside=512")
ax1.plot(gmap_cmb_256_k, gmap_cmb_256_k**2/(2*pi)**2 * gmap_cmb_256_pk, '-o', ms=2, c=c__256, label="nside=256")
ax1.set_ylim(min(np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10, 
             np.nanmin(gmap_cmb_1024_k[:20]**2/(2*pi)**2*gmap_cmb_1024_pk[:20])/10))
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor', title="Fourier power spectrum \nfrom the flat map $\\bf{generation}$ \nof the CMB power spectrum")

ax2 = plt.subplot(212)
ax2.loglog(gmap_cmb_1024_k, gmap_cmb_1024_k**2/(2*pi)**2 * gmap_cmb_1024_pk, '-o', ms=2, c=c_1024)
ax2.loglog(gmap_cmb_512_k, gmap_cmb_512_k**2/(2*pi)**2 * gmap_cmb_512_pk, '-o', ms=2, c=c__512)
ax2.loglog(gmap_cmb_256_k, gmap_cmb_256_k**2/(2*pi)**2 * gmap_cmb_256_pk, '-o', ms=2, c=c__256)
ax2.set_ylim(min(np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10, 
             np.nanmin(gmap_cmb_1024_k[:20]**2/(2*pi)**2*gmap_cmb_1024_pk[:20])/10))

power_axes_k(ax1, ax2)

fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(gmap_pwl_1024_k, gmap_pwl_1024_k**2/(2*pi)**2 * gmap_pwl_1024_pk, '-o', ms=2, c=c_1024, label="nside=1024")
ax1.plot(gmap_pwl_512_k, gmap_pwl_512_k**2/(2*pi)**2 * gmap_pwl_512_pk, '-o', ms=2, c=c__512, label="nside=512")
ax1.plot(gmap_pwl_256_k, gmap_pwl_256_k**2/(2*pi)**2 * gmap_pwl_256_pk, '-o', ms=2, c=c__256, label="nside=256")
ax1.set_ylim(min(np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10, 
             np.nanmin(gmap_pwl_1024_k[:20]**2/(2*pi)**2*gmap_pwl_1024_pk[:20])/10))
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor', title="Fourier power spectrum \nfrom the flat map $\\bf{generation}$ \nof the custom power spectrum")

ax2 = plt.subplot(212)
ax2.loglog(gmap_pwl_1024_k, gmap_pwl_1024_k**2/(2*pi)**2 * gmap_pwl_1024_pk, '-o', ms=2, c=c_1024)
ax2.loglog(gmap_pwl_512_k, gmap_pwl_512_k**2/(2*pi)**2 * gmap_pwl_512_pk, '-o', ms=2, c=c__512)
ax2.loglog(gmap_pwl_256_k, gmap_pwl_256_k**2/(2*pi)**2 * gmap_pwl_256_pk, '-o', ms=2, c=c__256)
ax2.set_ylim(min(np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10, 
             np.nanmin(gmap_pwl_1024_k[:20]**2/(2*pi)**2*gmap_pwl_1024_pk[:20])/10))

power_axes_k(ax1, ax2)

 

Compare flat projection from full-sky with flat map generation

Maps

vmax = max(abs(fmap_cmb_1024[0].min()), abs(fmap_cmb_1024[0].max()))
plt.title("flat projection from full-sky CMB map\n %.f x %.f deg$^2$ \n %d x %d pixels" 
          % (imgsize_arcmin/60, imgsize_arcmin/60, pixsize_1024, pixsize_1024))
plt.imshow(fmap_cmb_1024[0], cmap=cmap, vmin=-vmax, vmax=vmax)
plt.colorbar()

plt.figure()
vmax = max(abs(gmap_cmb_1024.min()), abs(gmap_cmb_1024.max()))
plt.title("flat map generation for CMB\n %.f x %.f deg$^2$ \n %d x %d pixels" 
          % (imgsize_arcmin/60, imgsize_arcmin/60, pixsize_1024, pixsize_1024))
plt.imshow(gmap_cmb_1024, cmap=cmap, vmin=-vmax, vmax=vmax)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fa20c10eb80>

vmax = max(abs(fmap_pwl_1024[0].min()), abs(fmap_pwl_1024[0].max()))
plt.title("flat projection from full-sky custom power-law map\n %.f x %.f deg$^2$ \n %d x %d pixels" 
          % (imgsize_arcmin/60, imgsize_arcmin/60, pixsize_1024, pixsize_1024))
plt.imshow(fmap_pwl_1024[0], cmap=cmap, vmin=-vmax, vmax=vmax)
plt.colorbar()

plt.figure()
vmax = max(abs(gmap_pwl_1024.min()), abs(gmap_pwl_1024.max()))
plt.title("flat map generation for custom power-law\n %.f x %.f deg$^2$ \n %d x %d pixels" 
          % (imgsize_arcmin/60, imgsize_arcmin/60, pixsize_1024, pixsize_1024))
plt.imshow(gmap_pwl_1024, cmap=cmap, vmin=-vmax, vmax=vmax)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fa20c4861c0>

Fourier power spectra

fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(fmap_cmb_1024_k, fmap_cmb_1024_k**2/(2*pi)**2 * fmap_cmb_1024_pk, '-o', ms=2, c=c_proj, label="projected full-sky map")
ax1.plot(gmap_cmb_1024_k, gmap_cmb_1024_k**2/(2*pi)**2 * gmap_cmb_1024_pk, '-o', ms=2, c=c_fgen, label="flat map generation")
ax1.set_ylim(min(np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10, 
             np.nanmin(gmap_cmb_1024_k[:20]**2/(2*pi)**2*gmap_cmb_1024_pk[:20])/10))
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor', title="CMB power spectrum")

ax2 = plt.subplot(212)
ax2.loglog(fmap_cmb_1024_k, fmap_cmb_1024_k**2/(2*pi)**2 * fmap_cmb_1024_pk, '-o', ms=2, c=c_proj)
ax2.loglog(gmap_cmb_1024_k, gmap_cmb_1024_k**2/(2*pi)**2 * gmap_cmb_1024_pk, '-o', ms=2, c=c_fgen)
ax2.set_ylim(min(np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10, 
             np.nanmin(gmap_cmb_1024_k[:20]**2/(2*pi)**2*gmap_cmb_1024_pk[:20])/10))

power_axes_k(ax1, ax2)

fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(fmap_pwl_1024_k, fmap_pwl_1024_k**2/(2*pi)**2 * fmap_pwl_1024_pk, '-o', ms=2, c=c_proj, label="projected full-sky map")
ax1.plot(gmap_pwl_1024_k, gmap_pwl_1024_k**2/(2*pi)**2 * gmap_pwl_1024_pk, '-o', ms=2, c=c_fgen, label="flat map generation")
ax1.set_ylim(min(np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10, 
             np.nanmin(gmap_pwl_1024_k[:20]**2/(2*pi)**2*gmap_pwl_1024_pk[:20])/10))
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor', title="custom power-law power spectrum")

ax2 = plt.subplot(212)
ax2.loglog(fmap_pwl_1024_k, fmap_pwl_1024_k**2/(2*pi)**2 * fmap_pwl_1024_pk, '-o', ms=2, c=c_proj)
ax2.loglog(gmap_pwl_1024_k, gmap_pwl_1024_k**2/(2*pi)**2 * gmap_pwl_1024_pk, '-o', ms=2, c=c_fgen)
ax2.set_ylim(min(np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10, 
             np.nanmin(gmap_pwl_1024_k[:20]**2/(2*pi)**2*gmap_pwl_1024_pk[:20])/10))

power_axes_k(ax1, ax2)

ax1 = plt.subplot(111)
ax1.loglog(fmap_cmb_1024_k, fmap_cmb_1024_pk, label="projected full-sky map")
ax1.loglog(gmap_cmb_1024_k, gmap_cmb_1024_pk, label="flat map generation")
ax1.set_xlabel("$k$")
ax1.set_ylabel("$P(k)$")
ax1.legend(handlelength=0, labelcolor='linecolor', title="CMB power spectrum");

ax1 = plt.subplot(111)
ax1.loglog(fmap_pwl_1024_k, fmap_pwl_1024_pk, label="projected full-sky map")
ax1.loglog(gmap_pwl_1024_k, gmap_pwl_1024_pk, label="flat map generation")
ax1.set_xlabel("$k$")
ax1.set_ylabel("$P(k)$")
ax1.legend(handlelength=0, labelcolor='linecolor', title="custom power-law power spectrum");